The characterization of toll‐like receptor repertoire in Pinna nobilis after mass mortality events suggests adaptive introgression

Abstract The fan mussel Pinna nobilis is currently on the brink of extinction due to a multifactorial disease mainly caused to the highly pathogenic parasite Haplosporidium pinnae, meaning that the selection pressure outweighs the adaptive potential of the species. Hopefully, rare individuals have been observed somehow resistant to the parasite, stretching the need to identify the traits underlying this better fitness. Among the candidate to explore at first intention are fast‐evolving immune genes, of which toll‐like receptor (TLR). In this study, we examined the genetic diversity at 14 TLR loci across P. nobilis, Pinna rudis and P. nobilis × P. rudis hybrid genomes, collected at four physically distant regions, that were found to be either resistant or sensitive to the parasite H. pinnae. We report a high genetic diversity, mainly observed at cell surface TLRs compared with that of endosomal TLRs. However, the endosomal TLR‐7 exhibited unexpected level of diversity and haplotype phylogeny. The lack of population structure, associated with a high genetic diversity and elevated dN/dS ratio, was interpreted as balancing selection, though both directional and purifying selection were detected. Interestingly, roughly 40% of the P. nobilis identified as resistant to H. pinnae were introgressed with P. rudis TLR. Specifically, they all carried a TLR‐7 of P. rudis origin, whereas sensitive P. nobilis were not introgressed, at least at TLR loci. Small contributions of TLR‐6 and TLR‐4 single‐nucleotide polymorphisms to the clustering of resistant and susceptible individuals could be detected, but their specific role in resistance remains highly speculative. This study provides new information on the diversity of TLR genes within the P. nobilis species after MME and additional insights into adaptation to H. pinnae that should contribute to the conservation of this Mediterranean endemic species.


| INTRODUC TI ON
Filtering-feeder marine invertebrates, such as bivalves, are exposed to a wide range of potentially pathogenic microbes including bacteria, viruses, and parasites. In most cases, host and pathogens have co-evolved over long evolutionary periods, allowing populations to cope effectively with the usual periodic epizootics (Altizer et al., 2003;Bean et al., 2013). This balance between pathogens and adapted host is currently challenged by ongoing anthropogenic changes, of which the introduction of new pathogens, global warming and pollution (Fey et al., 2015;Natalotto et al., 2015).
Newly introduced pathogens representing obvious threats for naive populations can cause dramatic demographic and genetic changes, eventually reducing biodiversity and altering ecosystem functioning and services (Elmqvist et al., 2003;Tompkins & Begon, 1999). From a biological conservation perspective, it is therefore necessary to quickly assess whether and to what extent a species is able to rapidly adapt to new infections and naturally recover (Mable, 2019). This knowledge is indeed crucial in determining the risk of global spread of pathogens, the magnitude of potential mass mortality events (MME) as well as the development of conservation strategies.
To this end, rapidly evolving immune genes are among the most relevant candidates genetic loci to study at first intention, in particular those encoding toll-like receptors (TLR) (Grueber et al., 2015;McCallum & Dobson, 2002;Turner et al., 2011). TLR indeed have a pivotal role in the recognition of a wide diversity of pathogens, through the binding of conserved pathogen-associated molecular patterns (PAMP) (Uematsu & Akira, 2008), as well as in the stimulation and regulation of both innate and adaptive immunity (Akira & Takeda, 2004;Barreiro et al., 2009). Previous research performed on different species has also shown that TLR genes have undergone pervasive or episodic positive selection (Fornarino et al., 2011;Grueber et al., 2015;Jiang et al., 2017;Tschirren et al., 2013), in that variations in the sequences of TLR genes could be associated with both pathogen resistance and disease tolerance (Chevillard et al., 2018;Destoumieux-Garzón et al., 2020;Imai et al., 2008;Jiang et al., 2017;Madan-Lala et al., 2017;Nahrendorf et al., 2021;Pila et al., 2016;Rakoff-Nahoum et al., 2004;Tschirren et al., 2013). Moreover, determining the diversity of adaptive loci may likely help to understand the evolutionary processes underlying adaptation to pathogens and provide biomarkers that could be used to identify eventually most adapted individuals (Mable, 2019;McCallum & Dobson, 2002). Different patterns of diversity could actually be expected according to the relative influence of selection, driven by one or multiple infectious agents for instance (Connallon & Clark, 2012;Hedrick, 1999), drift and demographic fluctuations. Selection can have antagonist effects by either maintaining or increasing diversity or reducing diversity by promoting the fixation of adaptive alleles in certain groups or local populations (Holderegger et al., 2006;Spurgin & Richardson, 2010). On the contrary, distinct scenarios of population size fluctuations can as well lead to contrasted measure of diversity (Sutton et al., 2015;Taylor & Jamieson, 2008).
For instance, excess of low-frequency variants would sign a strong long-lasting bottleneck, while an excess of intermediate-frequency variants would be related to a recent bottleneck of weak intensity, or to a balancing selection.
The diversity of TLRs has been mainly studied in vertebrates, globally highlighting potentially large number of haplotypes as well as high number of segregating sites . Noteworthy, in vertebrates, dN/dS ratios are generally less than unity, indicating that selection rather suppresses changes in protein sequences (i.e. negative or purifying selection) (Kryazhimskiy & Plotkin, 2008).
Regarding mollusc taxa, less knowledge is available, but a high level of polymorphism has been evidenced overall the transcriptome of the bivalve species Macoma balthica (Pante et al., 2012).
The possibility of obtaining a fairly complete repertoire of full-length coding DNA sequences is now possible thanks to nextgeneration sequencing, making it possible to understand how these genes are eventually shaped by selective pressures, of which pathogens.
The Mediterranean endemic pen shell Pinna nobilis (Linnaeus, 1758) is an iconic species living in Posidonia oceanica meadows, in the coastal fringes from 0 m to 60 m depth. Pinna nobilis plays a central ecological role, reducing the turbidity of the water as a filtering organism, sheltering a number of symbiotic species and providing a hard substrate for the development of many epibiotic species.
In 1992, P. nobilis was declared a protected species with the category of vulnerable (European Council Directive 92/43/EEC), mainly to counteract intensive captures and habitat loss. These conservation measures proved efficient in several regions, for instance in France, Croatia, Spain and Italy, where populations regained vitality (Marrocco et al., 2019;Siletic & Peharda, 2003;Trigos-Santos & Vicente, 2018). The most recent research on the dynamics and structure of the P. nobilis population indeed demonstrated high effective population sizes, overall genetic homogeneity between populations (i.e. based on F ST calculations) and a potentially high level of connectivity Peyran et al., 2021;Wesselmann et al., 2018).
Unfortunately, the species has been experiencing mass mortality events (MME) since 2016, attributed to a multifactorial disease mainly conveyed not only to the protozoan parasite Haplosporidium K E Y W O R D S adaptation, critically endangered, fan mussel, genetic diversity, Haplosporidium pinnae, hybrids, introgression, pathogens, toll-like receptors

T A X O N O M Y C L A S S I F I C A T I O N
Conservation genetics, Parasitology, Population genetics pinnae but also to the bacteria Vibrio mediterranei and Mycobacterium sp., among potentially yet uncharacterized other pathogens, such as viruses and micro-eukaryotes (Hamelin et al., 2019;Künili et al., 2021;Lattos et al., 2020;Scarpa et al., 2020). These MME led to dramatic demographic bottlenecks, throughout the species distribution area (Andree et al., 2020;Catanese et al., 2018;Lattos et al., 2020;Šarić et al., 2020;Vázquez-Luis et al., 2017), that are expected to critically reduce the genetic diversity, increase population fragmentation and probably negatively influence the natural rescue effect through larval connectivity, due to reduced reproduction success (Cowen et al., 2007). Therefore, although few populations have recently been reported to be safe in some places (Cinar et al., 2021;Peyran et al., 2022), the vast majority of populations have been affected by those pathogens. As a consequence, the species has been included in the IUCN Red List as Critically Endangered (Kersting et al., 2019) Green & Speck, 2018;Melillo et al., 2018;Wang et al., 2020) which would be a glimmer of hope for the conservation of the species.
In this study, our main objective was to bring new insights on the molecular mechanisms involved in the adaptation to H. pinnae epizootic (Salis et al., 2022). We chose a candidate approach focussing on TLR genes to test the hypothesis that TLRs are involved in such adaptation. Towards this aim, we obtained individuals (P. nobilis, P. rudis and their hybrids) from four distant geographical regions covering the largest part of all three taxa geographic ranges, sampled after epizootic diseases, in which H. pinnae was mainly involved. Therefore, we here described the genetic polymorphism and diversity at 14 predicted TLRs identified in the shotgun genome of P. nobilis (Bunet et al., 2021), which displayed classical extracellular leucine-rich repeats (LRRs) and intracellular Toll/Interleukin-1 receptor (TIR) domains (Gay & Gangloff, 2007). Then, we assessed whether and to what extent the observed genetic polymorphism could be associated with better MME survival.

| Pinna spp. samples
No experimentation involving removal, dislocation or killing of Pinna spp. individuals was performed. A nonlethal sampling method was specifically developed to avoid the manipulation of animals. A small piece (14 mm 2 , ~20 mg) of the mantle tissue was quickly excised using biopsy forceps before the fan mussel closes its valves. The piece of tissue was transferred in a tube filled with 95% ethanol and stored at −80°C until use.
Overall DNA or tissue samples of 43 adult individuals were obtained from Greece, Italy, Spain and France (Table 1, Table S1).
Among those, 38 were identified as Pinna nobilis, 1 as Pinna rudis and 4 as P. nobilis × P. rudis hybrids. These individuals were sampled during or after epizootic disease mainly due to Haplosporidium pinnae. Dead or moribund individuals at the time of sampling are termed 'sensitive' (S). For convenience, the others that were alive without obvious symptom of infection are termed 'resistant' (R) (see Section 4). The Pinna rudis and hybrids individuals did not present any signs of disease and are termed 'naturally resistant' (NR). Dead specimen presented heavy lesions in the digestive gland connective tissue and signs of degenerative process in epithelial tissue. Moribund specimen presented half-open valves and a retracted body ( Figure S1).
The presence of several pathogens has been investigated as described below, and using a diagnostic quantitative PCR as well for the detection of H. pinnae, as previously described (López-Sanmartín

| Read mapping and phasing
High-quality paired reads (Q > 40) with a read length higher than 150 nucleotides were aligned against the reference shotgun genome, as well as the partial sequences of Haplosporidium pinnae coding for the small subunit ribosomal RNA (LC338065), Mycobacterium sp. coding for hsp65 (AY379077) and ITS (MN637877), Vibrio sp.
coding for ATPaseA (MN201557), Perkinsus mediterraneus, P. atlanticus and P. marinus coding for the small subunit ribosomal RNA (AY486139, AF509333 and AF126013, respectively), using bwamem ). Sam files were processed using samtools ) in order to obtain bam files of aligned reads. The observed sequencing depth at coding DNA sequence was at least 10×.
Reads that aligned the sequences of 19 neutral microsatellites, 14 TLRs and the cytoplasmic dynein (i.e. used as a nonimmune control gene), were phased using the 'phase' module of samtools with default parameters, and haplotypes were then inferred using SPADES (parameters: -k 21,33,55,77,99,127 --careful --only-assembler) (Bankevich et al., 2012). All reconstructed sequences were visually checked using IGV (Thorvaldsdóttir et al., 2013). SNP present at less than 20% of the reads was not considered. Finally, each TLR haplotype was phased a second time using the phase module of DNAsp v6.0 allowing events of recombination.

| TLR structure, phylogenetic analyses and ligand-binding sites analyses
Once the haplotypes of predicted TLR gene-coding sequences were obtained, the coding DNA sequences and predicted protein sequences were deduced using Augustus (Stanke & Morgenstern, 2005) using default parameters. Then, the protein structures were crosschecked using SMART , PROSITE, TMHMM-2.0 (Krogh et al., 2001) and BLASTP in order to provide information about the presence of leucine-rich repeat motifs (LRR) and Toll/IL-1 resistance (TIR) domains which are functionally important and transmembrane domains.
Then, sequence consensus for each TLR protein was obtained from multisequence alignment using EMBOSS Cons available online (https:// www.ebi.ac.uk/Tools/ msa/emboss_cons/), for both P. nobilis and P. rudis. Consensus sequences were used to model the 3D structure, using I-TASSER (Yang et al., 2015) and PyMOL (http://www.pymol. org/pymol), and to draw a phylogenetic tree depicting the relationships between TLR proteins based on multiple sequence alignments, using MEGA 11 and visualized using iTOLv6 ).

| Genetic diversity and structure
Neutral microsatellites and TLR haplotypes were aligned with MEGA11, using CLUSTAL and MUSCLE (codon), respectively. Singletons of TLR sequences were removed. These alignments were used in DNAsp (Rozas et al., 2017) to determine haplotype diversity, which were translated into genotype data. Regarding TLR haplotypes specifically, the following diversity statistics were assessed using DNAsp (Rozas et al., 2017): the number of polymorphic sites (S), the average number of nucleotide differences among all alleles (k), the nucleotide diversity (π), the number of nonsynonymous (dN) and synonymous (dS) changes. Prior to the genetic structure analysis, ARLEQUIN v3.5.1.2 (Excoffier & Lischer, 2010) was also used to check that the retained loci were not in linkage disequilibrium, to assess the distribution of genetic diversity by an analysis of molecular variance (AMOVA) and to calculate the pairwise genetic differentiations (F ST ) among each population. The corrected significance threshold for multiple tests was set using the Benjamini-Hochberg correction procedure (Benjamini & Hochberg, 1995). Subsequently, genetic structures were depicted using principal component analysis (PCA) using R. were implemented with a total of 150,000 iterations, a burn-in period of 50,000 and 20 pilot runs of 10,000 iterations each and did not yield any significant results. The analysis using FUBAR (Murrell et al., 2013) was performed once potential recombination events were assessed for each TLR locus using GARD algorithm (Kosakovsky Pond et al., 2006), also accessible from the DataMonkey server. This is a recommended preprocessing step for selection inference, so that any recombination effects could be considered in the FUBAR model.

| Assessment of selection acting at TLR loci
Only sites under putative selection with a posterior probability higher than 0.9 were retained for further analysis. Finally, DAPC analysis was performed using R, in order to depict the relative contributions of SNP retaining the three to four clusters of individuals (i.e. according to their phenotypes), and with n.pca = 6; n.da = 4.

| Pinna spp. samples and infectious status
The 43  In both P. nobilis and P. rudis species, we found one or several members of TLR-1, -2, -3, -4, -6, -7, -13 and toll-like families (Table S2). The phylogenetic relationships between TLR coding sequences of both species, as well as the protein domains of each TLR, are depicted in Figure S2.
Across the 14 P. nobilis TLR loci, a total of 228 haplotypes were observed with a mean haplotype per loci of 16.4 ± SD 8.1 and a mean haplotype diversity (H d ) of 0.771 ± SD 0.166. Overall, one or a few haplotypes predominate at each TLR locus (Figure 3, Figure S3). The density of segregating sites is rather homogenous across loci, thus between both the extra and intracellular domains of the protein, with a mean of roughly 8.2 ± SD 5.8 substitutions per 1000 nt, except for the TLR-7 in which 27 substitutions per 1000 nt was observed ( Table 2).
The dN/dS ratio observed at the full-length sequences ranged between 0.5 and 3.8 (mean = 1.4 ± SD 1.0) and was mostly higher at the extracellular domain (mean = 2.1 ± SD 2.1) compared with that of the intracellular domain (mean = 0.6 ± SD 0.6), which reveals a greater polymorphism in the protein region harbouring binding sites.
The observed mean number of translated proteins was 11 ± SD 6.2 (Table 2 and Figure 1).
The polymorphism observed at TLR loci was found like that observed for the nonimmune comparative control gene (i.e. 15 haplotypes, H d = 0.714 and S/1000 nt = 9.4). However, for dynein, the dN/ dS ratio was 0.2, that is three and seven times lower than the dN/ dS ratios observed at intracellular domain and at the full-length sequence of TLR, respectively.
As expected for TLR genes, the genetic polymorphism observed between orthologue sequences was rather high, with a mean density of segregating sites of 42.1 ± SD 12.6 sites per 1000 nt (Table 2).

| Basic genetic diversity and distribution
Similar levels of genetic and gene diversities have been detected among populations as well as among phenotype groups, for both microsatellites and TLR loci. Overall, no sign of departure from Hardy-Weinberg equilibrium was observed since F IS were not significant, and specifically no excess of heterozygotes could be observed (  (Figure 2, Figure S3).

| Assessment of selection at TLR loci
3.4.1 | Tests of departure from neutrality As no genetic structure could be obviously evidenced, and to prevent from bias resulting from low population size, we considered the whole sequences as coming from one panmictic population.

| Pattern of nonsynonymous versus synonymous substitutions
Considering full-length sequences, 7 TLR presented a dN/dS ratio higher than 1, of which three comprised between 1.5 and 3.

| Median-joining networks of haplotypes
Overall, we observed clear star-like phylogenies, with few highfrequency nodes from which derive eventually rare haplotypes, which could be consistent with a population expansion model and/ or a recent positive selection ( Figure S3). This pattern of population expansion may reflect the increase of individuals, which is not necessarily accompanied by recovery in the loss of genetic  Figure S3).

| Phylogeny-based test of selection and discriminant analysis of principal component (DAPC)
Before the analysis using FUBAR, we assessed the putative recombination sites within TLR sequences, using GARD. A total of eight recombination sites were detected across six TLR sequences, mostly consisting in one recombination site detected in four TLR sequences and two recombination sites in TLR-7 and a protein toll-like (data not shown).
Among the 1481 polymorph sites observed across TLR loci (considering both P. nobilis and P. rudis), FUBAR identified 395 sites (26.6%) putatively under selection, with a posterior probability of .9.
Most of these (95%) were negatively selected sites, and most of the positively selected sites (47 out of 53 sites) were found at the extracellular domain of the proteins (Figure 1).
Complementary to FUBAR analysis, DAPC performed using either all or the putatively selected SNPs evidenced that the group of naturally resistant phenotype was obviously separated from the two others, with a probability of assignment of individuals to this cluster of 100% (data not shown). On the contrary, the groups of resistant and sensitive phenotypes overlapped, though few resistant individuals were remotely distributed from the others (Figure 4).
When the entire set of SNPs was used in the DAPC, almost all SNP of TLR-4 (i.e. Contig 17,440), one SNP of TLR-6 (Contig 48,600, snp771) and one SNP of TLR-1 (Contig 21,812, snp1030) mainly contributed to the cluster according to axis 1, which separate the group of naturally resistant individuals from other groups.

| DISCUSS ION
Determination of the evolutive capacity of threatened species is among the major keystone actions towards their conservation. Among the genes to be monitored at first intention are immune genes which act as sentinels against infectious agents. In the context of high P. nobilis mortality, our aim was to describe the current genetic polymorphism of the rapidly evolving toll- F I G U R E 5 TLR-7 polymorphism and 3D structure. (a) Sequence alignments of several regions of TLR-7 of P. rudis origin within introgressed P. nobilis. Consensus sequences of P. nobilis and P. rudis have been included for comparison. PSS, amino acids under positive selection as detected by FUBAR. (b-d) the 3D models have been obtained using I-TASSER software from P. nobilis and P. rudis TLR-7 consensus sequences, represented in green and blue, respectively. (b) superimposed TLR-7 protein at site 1 with 06S ligand; (c) P. nobilis and P. rudis TLR-7 at site 2 with XG1 and RX8 best-fit ligands, respectively; (d) superimposed TLR-7 protein around site 733. Divergent amino acid substitutions are represented in pink. Greece, up to a year for those of Italy), but finally died, which may be attributed to other pathogens or factors as well. During MMEs, no regular monitoring of disease progression at the individual level could be done, so it was not possible to clearly compare the level of disease at the specimen level, nor the interindividual stress response over the course of infection. Given our samples and the lack of standardized sampling protocol, such an assessment could not be made neither. As a consequence, strictly speaking, these P. nobilis cannot be obviously defined as either resistant or tolerant, precisely to H. pinnae infection, despite the eventual inability to survive. In our analysis, for convenience, these P. nobilis were considered resistant.
Although we cannot exclude some assignment errors in the resistant group, we assumed that field observations were nevertheless good indicators of signs of resistance. At the molecular level, we were able to demonstrate an obvious difference between the sensitive and resistant groups since introgression at the TLR-7 locus was observed only in the group of resistant P. nobilis. Additionally, three individuals of P. nobilis carrying the original full-length TLR-7 of P. rudis appeared as naturally resistant as the hybrids and P. rudis (i.e. no infection detected), suggesting that TLR-7 was indeed probably involved in disease resistance.
This result was obtained considering a rather moderate sample size from a field study. Actually, the lack of surviving populations is a critical issue throughout the species distribution, discommoding the collection of larger number sizes of individuals, a situation that has recently been discussed (Salis et al., 2022). Nonetheless, despite the small number of TLR-7-introgressed naturally resistant P. nobilis of our dataset, this result appears highly consistent with the recent research evidencing that TLR-7 can efficiently hamper the development of intracellular protozoan parasites (Baccarella et al., 2013;Fouzder et al., 2019;Heni et al., 2020;Regli et al., 2020). Therefore, we argue that it could rapidly help to develop conservation strategies for the species survival.
To clarify the relative importance of TLR-7 introgression in the resistance to H. pinnae, further laboratory-controlled studies should be performed. On the contrary, introgression should rapidly be more investigated in future monitoring studies as well, focussing on the TLR-7 genetic composition of Pinna individuals collected either from mortality events or surviving populations.
This knowledge will inform about the evolutionary potential of introgressed P. nobilis and the contribution of these individuals to future generations. Probably, most of the individuals identified as resistant had the same null contribution to the susceptible ones (Hasik & Siepielski, 2022). Whether the introgressed specimen contributes more to the next generation than the susceptible and the other resistant individuals, in the context of MME, needs to be assessed in future monitoring studies.

| TLR validation and subtypes identification
Except for TLR-4 (i.e. contig 473) and TLR-4/2/13 (i.e. contig 17,440), all coding DNA sequences of our contigs consisted in one large exon, a feature that has been reported in other molluscs (Juhász & Lawton, 2022). All TLR coding sequences likely encoded functional receptors, as we did not observe any stop codon, and since the expected LRR and TIR motifs were always detected. The protein annotation performed using protein BLAST evidenced that the TLR repertoire we studied was composed of cell surface TLRs (i.e. TLR-1/ toll-like protein, TLR-2, TLR-4 and TLR-6) and endosomal TLRs (i.e. TLR-3, TLR-7/Tollo and TLR-13) and some levels of gene redundancy. This is consistent with recent evidence that immune genes are often tandemly duplicated in molluscs and specially in bivalves (Batista et al., 2019;Peng et al., 2020;Zhang et al., 2013). Such an observation was, as well, observed on one contig (i.e. contig 7594) carrying the related TLR-13 and TLR-13/3.
Finally, the 14 TLR loci presented here are part, yet probably nearly complete, of the TLR repertoire of the two species. Indeed, 21 contigs were annotated as TLR in the shotgun genome of which five encoding partial TLR sequences (i.e. either LRR or TIR for instance) and two that aligned a high proportion of the same reads, which would have strongly biased the analysis.

| Variation of TLR genes
Haplotype networks and the negative values of Tajima's D overall indicated that the diversity, at most TLR loci, was made of few highfrequency haplotypes and a few to high number of rare alleles. The number of haplotypes and the level of polymorphism were similar to that reported in several vertebrate and invertebrate species . Densities of segregating sites observed across TLR genes were similar to that of the nonimmune control gene and are consistent with previous results obtained from the transcriptome sequencing in the bivalve Macoma balthica (Pante et al., 2012), which would indicate that we likely did not overestimate the level of polymorphism.
Noteworthy, TLR loci were not equally polymorphic. The observed differences in polymorphism were overall consistent with expectations in that cell surface TLRs, which bind an expanded set of ligands of different nature, exhibit high allele diversity, compared with endosomal TLRs which sense rather conserved nucleic acid structures (Georgel et al., 2009). In our study, we observed a rather unexpected high number of haplotypes at the vesicular TLR-7 which usually presents a low diversity, at least in vertebrates. This result could be consistent with a relaxed selection acting at this locus, reflecting the high diversity of intracellular pathogens that could possibly be encountered in the environment (Künili et al., 2021;Lattos et al., 2020). Another nonexclusive explanation would be that introgression brings high levels of diversity at that locus (Fraïsse et al., 2016), which could be transmitted all in one, or partially to successive generations due to intragenic recombination. Recombination would then be a driver of diversity leading to the production of a subset of haplotypes, which then might evolve differentially than the main haplotypes. Such a contribution of recombination has already been evidenced in immune genes, including TLR (Gösser et al., 2019;Oren et al., 2019;Schaschl et al., 2006). To our knowledge, the diversity of TLR-7 has not been described in bivalves, and a high level of TLR-7 has just been observed in chicken (Bulumulla et al., 2011;Świderská et al., 2018). Conversely, the two other types of vesicular TLRs (i.e. TLR-13 and TLR-3) sense more structurally conserved double-stranded bacterial RNAs.

| Assessment of selection towards TLR loci
No population structure of P. nobilis could be detected, neither using neutral microsatellites nor TLR loci, and high and similar levels of genetic diversity were also revealed for both markers, consistent with previously published research (Peyran et al., 2021;Wesselmann et al., 2018). Moreover, considering the overall haplotypes topologies, as well as the negative Tajima's D, our results are in agreement with the hypothesis of a panmictic population that has been experiencing a recent postbottleneck expansion (Sanna et al., 2014) and suggest that demographic processes overweigh, or hamper the detection of a small effect of selection associated with the resistance to H. pinnae.
The excess of nonsynonymous mutations and high heterozygosity rates detected is, as well, congruent with at least a transient balancing selection, acting at some of the TLR loci studied Minias & Vinkler, 2022). In fact, observing a clear correlation between genetic traits and a specific pathogen resistance is difficult when the populations being tested likely have contrasting life histories with respect to the microbial communities whose pathogens they have met. Accordingly, pathogens are expected to maintain high genetic polymorphisms in hosts' populations, again interpreted as balancing selection (Llaurens et al., 2017;Quéméré et al., 2018). The proposed mechanism for maintaining such a high diversity is based on either the advantage of heterozygotes (i.e. overdominance), the advantage of rare alleles (Browne & Karubian, 2018;Stefan et al., 2019), fluctuating selection (Bell, 2010;Quéméré et al., 2018)  In this study, we also report that hybrids are indeed significantly introgressed and that more than 18% of the P. nobilis of the dataset originated from ancient hybrids. Noteworthy, all P. nobilis we found introgressed at TLR loci were identified as resistant and carried a TLR-7 allele of P. rudis origin. As discussed above, TLR-7 is increasingly recognized as a potent modulator of protozoan infection, capable of impeding parasite growth (for an extensive review see Fouzder et al., 2019). Therefore, the fact that TLR-7 is conserved in MMEsurviving P. nobilis could be interpreted as an obvious a sign of adaptive introgression (Edelman & Mallet, 2021;Gouy & Excoffier, 2020;Hedrick, 2013), but once again, there is a strong indication that the full-length TLR-7 of P. rudis is required for the effective resistance to H. pinnae.
Once activated, TLR-7 homodimers activate the interferon regulatory factors (IRF) signalling pathway and trigger both autophagy and the production of type I interferon and other inflammatory cytokines (Petes et al., 2017).
Structurally, TLR-7 possesses two binding sites for a guanosine at site 1 and a single-stranded RNA (ssRNA) at site 2. Both binding sites are situated at either side of a Z-loop region which is also involved in the ssRNA recognition. Recent research evidenced that the first site was essential for TLR-7 activation, while the second one only contributes to ssRNA-induced activation. The TLR-7 is thus synergistically activated in response to guanosine and ssRNA (Zhang et al., 2016(Zhang et al., , 2018. This feature of TLR-7 likely explains why the 4 P. nobilis introgressed with a partial TLR-7 of P. rudis could not effectively resist the infection with H. pinnae, whereas those with a full-length P. rudis TLR-7 were not affected. Regarding TLR-6 and TLR-4, whether these TLR are indeed more or less involved in some adaptation to H. pinnae remains highly speculative and should require specific studies. Nonetheless, these TLR, once bound to their ligands, activate MyD88 and NFkB signalling pathways, hence triggering a pro-inflammatory response (Krishnegowda et al., 2005;Molteni et al., 2016) which may potentialize the specific TLR-7 response.
TLR-6 typically binds diacylated lipopeptides (Drage et al., 2009) and glycoproteins glycophosphatidylinositol (Takeuchi et al., 2001) and has been implicated in the response and potential resistance to several pathogens, from viruses (Chen et al., 2015) to bacteria (Marinho et al., 2013) and protozoan (Leoratti et al., 2008;Netea et al., 2008). Functionally, TLR-6 dimerizes with TLR-2, directing specific interaction with a pathogen molecule, and hence triggers essentially an IL-6 and TNFα mediated response (Krishnegowda et al., 2005). Interestingly, several research highlighted the role of TLR-6 in mitigating Mycobacterium infection (Marinho et al., 2013) to which P. nobilis is regularly exposed to. Certain adaptive changes towards the bacteria could confer a better ability to cope with the H. pinnae infection as well.
TLR-4 is mostly a sensor of the bacterial lipopolysaccharide (LPS) (Qureshi et al., 1999), although a wide range of bacterial, viral and protozoan ligands have also been identified (Molteni et al., 2016;Uematsu & Akira, 2008), indicating that TLR-4 participates in the response to a wide range of pathogens.

| CON CLUS ION
This preliminary research supports the role of TLR-7 in the resistance of P. nobilis to H. pinnae and the possible influence of TLR-6 and TLR-4, although of much smaller effect. Further field and laboratory studies are needed to complete and clarify the role of these TLRs in disease resistance and tolerance. In this aim, a systematic and extensive characterization of health and infectious status will be required. From field research, it would be necessary to sample more P. nobilis that survived at least one H. pinnae-driven MME and assess whether they are indeed introgressed by full-length P. rudis TLR-7. Similarly, systematically determining young recruit's haplotypes at TLR-7 locus could, as well, help predict how resilient a population would be in the face of an H. pinnae epizootic. At the laboratory, controlled experimental infections, using genetically well-characterized individuals would bring the evidence that TLR-7 is indeed involved in the adaptation to H. pinnae and could also allow to assess the level of associated stress response (Box et al., 2020;Lattos et al., 2023;Natalotto et al., 2015). Moreover, prophylactic strategies to prevent infection or enhance appropriate physiological responses should be tested. Current research on this topic seems promising, by challenging the innate immune system of the host with TLR-7 agonists which could modulate the innate immune system, thus either reducing tissue inflammation, help eliminate parasites, or prevent entry of the parasite in the host (El Hajj et al., 2018;Hamie et al., 2021). However, such experiment would probably be difficult to be performed and expensive. Finally, the fact remains that to date, no parasitic DNA has been demonstrated in P. rudis, hybrids or the P. nobilis introgressed with P. rudis full-length TLR-7, evidencing the likely pivotal role of this protein. Whether other pattern recognition receptors, eventually of P. rudis origin, are also involved in the protection against H. pinnae infection remains to be determined. The ongoing molecular research will probably bring additional insights in the mechanism of resistance.

ACK N OWLED G M ENTS
We thank Dr. Santiago V. Jiménez-Gutiérrez for providing us with some tissue samples of hybrids from Spain. We also thank an anonymous reviewer for helpfull comments.